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Abstract. 


The  determination  of  source  signature  is  a  major  calibration  problem  in  reflection  seismology . 
This  “deconvolution”  problem  is  conventionally  approached  by  way  of  statistical  methods,  by 
direct  measurement,  or  by  the  location  of  a  clean  reflection  in  an  otherwise  quiet  part  of  a 
reflection  section.  We  show  that  a  quasi-impulsive,  isotropic  point  source  may  be  recovered  simul¬ 
taneously  with  the  velocity  profile  from  reflection  data  over  a  layered  fluid,  in  linear  (perturbation) 
approximation.  Our  approach  is  completely  deterministic,  and  does  not  depend  on  the  presence  of 
an  isolated  reflection  in  a  quiet  part  of  the  section,  as  we  illustrate  with  a  numerical  example. 
After  describing  the  algorithm  and  a  numerical  implementation,  we  give  a  complete  mathematical 
treatment,  which  shows  that  our  estimates  of  source  wavelet  and  velocity  profile  are  stable  in  a 
certain  sense.  Because  of  this  stability  property  we  conjecture  that  our  approach  to  simultaneous 
estimation  of  source  and  medium  parameters  actually  applies  to  a  much  broader  class  of  models 


than  that  treated  here. 
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1.  Introduction. 

A  major  calibration  problem  in  reflection  seismology  is  that  of  source-signature  identification. 
The  seismic  section  depends  not  just  on  the  mechanical  structure  of  the  earth,  the  elucidation  of 
which  is  the  aim  of  reflection  seismology,  but  also  on  the  time  dependence  and  spatial  distribution 
of  the  seismic  energy  source.  The  separation  of  these  two  factors  (earth  structure,  energy  source)  is 
often  regarded  as  a  deconvolution  problem  (“source-signature  deconvolution”)  and  is  commonly 
attacked  with  statistical  tools  (predictive  deconvolution,  ARMA  modeling,  maximum  likelihood 
estimation),  which  aim  simultaneously  to  suppress  other  sorts  of  seismic  “noise”  (notably  multiple 
reflections).  These  widely-used  methods  have  been  criticized  as  being  based  on  unwarranted 
assumptions  about  both  source  and  earth  structure  (Ziolkowski,  1984).  An  alternate  method  of 
source  calibration  is  actual  measurement  of  the  direct  wave,  or  use  of  a  strong  reflection  in  an  oth¬ 
erwise  relatively  quiet  part  of  the  section. 

In  this  paper  we  investigate  the  possibility  that  the  time-dependence  of  an  isotropic  point- 
source  might  be  determined  directly  (and  deterministically)  from  the  section,  even  when  no  isolated 
reflection  event  can  be  located.  In  effect,  we  ask  whether  the  source  parameters  can  be  determined 
simultaneously  with  the  earth  model.  We  show  that,  under  some  restriction^-the  answer  is  “yes.” 

This  possibility  of  codetermination  of  source  and  velocity  arises  from  the  different  dependence 
of  their  data  influence  on  offset:  that  is,  the  effects  of  changing  the  source  wavelet,  respectively  the 
velocity  structure,  move  out  and  scale  differently.  These  effects  may  be  separated  algebraically; 
there  results  an  equation  linking  the  data  section  with  the  velocity,  sampled  at  several  different, 
but  related,  depths.  The  wavelet  has  been  eliminated  from  this  relation  altogether.  Once  the  velo¬ 
city  has  been  determined,  the  wavelet  may  be  extracted  easily.  On  the  other  hand,  functional- 
delay  equations,  such  as  the  one  we  derive  here  for  the  velocity,  are  not  commonplace  in  the 
mathematics  of  seismology,  and  it  is  not  entirely  obvious  that  such  equations  can  be  solved  in  a 
reasonable  sense.  Accordingly  we  give  a  complete  theoretical  treatment,  ensuring  that  solutions 
can  be  found  computationally  as  well.  We  implement  an  algorithm  based  on  our  theoretical 
developments,  and  show  by  numerical  illustration  that  wavelet  and  velocity  may  indeed  be 
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separately  determined  even  when  their  influence  is  coexistent  in  time. 

To  derive  this  result,  we  have  imposed  some  hypotheses  on  both  the  earth  model  and  the 
energy  source.  Those  concerning  the  earth  model  are: 

(Ml)  that  the  earth  structure  varies  only  with  depth  (layered  medium); 

(M2)  that  the  earth  is  a  linearly  elastic  fluid,  with  only  the  sound  velocity  varying  with 
position  (depth); 

(M3)  that  the  variation  of  sound  velocity  with  depth  is  somewhat  smooth. 

Concerning  the  source,  we  have  assumed  that 

(51)  the  spatial  distribution  of  the  source  is  point-like,  to  adequate  approximation; 

(52)  the  source  radiation  pattern  is  isotropic,  so  that  the  source  is  described  by  a  single 
function  of  time  (wavelet); 

(53)  the  source  wavelet  is  quasi-impulsive,  i.e.  differs  from  the  Dirac  delta  function  by  a 
square-integrable  function. 

None  of  these  assumptions  (with  the  possible  exception  of  (Si))  is  valid  for  an  accurate  model 
of  the  typical  seismic  reflection  experiment:  The  earth  is  non-layered  and  inhomogeneous  on 
interesting  length  scales  and  supports  shear  (as  well  as  compressional)  waves.  Also  the  seismic 
source  often  has  an  anisotropic  radiation  pattern,  and  (most  important)  is  bandlimited  (not  close 
to  delta). 

It  seems  clear  that  most  of  these  assumptions  could  be  relaxed,  at  least  to  some  extent.  The 
most  important  and  involved  step  is  the  introduction  of  non-impulsive  sources.  It  is  plausible  on 
the  basis  of  velocity  analysis  (and  is  carefully  justified  on  theoretical  grounds  in  Santosa  and 
Symes  (1986))  that  bandlimited  point-source  data  (with  known  source)  determine  a  layered  velocity 
structure,  provided  that  the  target  structure  is  sufficiently  rich  in  reflectors  (thus  is  sufficiently 
non-smooth).  We  expect  to  be  able  to  combine  these  bandlimited  inversion  ideas  with  those 
presented  in  this  paper  to  codetermine  bandlimited  sources  as  well.  We  will  further  discuss  the 
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extent  to  which  the  restrictions  M1-M3  and  S1-S3  might  be  relaxed  in  the  concluding  section. 

We  note  that  Canadas  and  Kolb  (1986)  give  numerical  evidence  that  the  velocity  and  source 
wavelet  (and,  to  some  extent,  the  density)  of  a  layered  fluid  may  be  recovered  from  simulated 
seismic  reflection  data.  They  do  not,  however,  give  a  theoretical  justification  for  their  procedure. 
Ramm  (1987)  has  given  a  theoretical  result  for  quasi-impulsive  sources,  which  differs  from  ours  in 
that  it  makes  use  of  the  low-frequency  (rather  than  high-frequency)  asymptotics  of  the  response, 
hence  is  intrinsically  restricted  to  the  quasi-impulsive  case. 

The  book  by  Lavrent’ev  et  al.  (1986)  has  recently  come  to  our  attention  (January  1988).  In 
Section  2.2  the  authors  consider  a  problem  similar  to  that  studied  here.  Their  work  is  similar  to 
ours  in  that  they  also  find  that  the  recovery  of  first-order  perturbations  in  medium  and  source 
parameters  hinges  on  the  solution  of  a  multiplicative-delay  equation  like  (3.2)  below.  It  differs 
from  ours  in  that  they  do  not  use  the  plane- wave  decomposition,  and  restrict  their  attention  to 
perturbations  about  a  homogeneous  background  model. 


The  model  described  above  is  quantified  in  the  following  boundary  value  problem,  connecting 


p(*) 

-  fluid  density 

X(x) 

-  bulk  modulus 

«(z,  0 

-  particle  displacement  field 

/(<) 

-  source  wavelet: 

P(x )  0  =  V(X(z)V-  «(x,  t)) 

in  {(z ,  t ):  z3  >  0} 

X(z)  V  •  u(x,  t)  =  -  f{t)d(x1,x2) 
for  z3  =  0 

u(z,  t)  =  0.  t  «  0 


We  write  z  =  z3,  and  note  that  (Ml)  is  stated  p  =  p{z),  X  =  X(z).  We  assume  that  X  =  const.  , 

so  that  e(z)  =  X2p(z)  2  parameterizes  the  medium  (M2). 

We  remark  that  this  is  not  a  significant  restriction:  the  problem  of  simultaneous  determina¬ 
tion  of  p  and  X  from  known  (and  impulsive)  source  data  is  well  understood  (see  e.g.  Santosa  and 
Symes  (1987)  and  other  work  cited  there).  Extension  of  our  present  results  to  this  more  general 
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case  is  routine.  We  have  chosen  to  present  the  result  for  \=const.  merely  to  avoid  obscuring  the 
novelty  of  the  present  paper  with  irrelevant  complexities. 

The  hypothesis  (S3)  is  that 

/(0«*(0  +  /i(0 

with  / 1  square-integrable.  Causality  is  expressed  by  the  requirement  that  f  i(t)  =  0,  t  <  0  and  we 
also  assume  that  fi(t)  =  0,  t  >  fm„.  Consequently,  /  defines  an  invertible  convolution  operator 
(on  L2[ 0,  T]  for  any  T  >  0). 

We  idealize  the  seismogram  as  the  surface  trace  of  the  vertical  displacement  field 
u  —  (vJt  uf  ut)  viewed  as  a  function  of  the  velocity  profile  c  and  the  wavelet  /  : 

*(«»/]  =  «*(*  =  0)  • 

We  shall  use  the  commonplace  notation  of  5(  • )  to  denote  perturbations  of  ( ■ );  the  distinction  from 
the  Dirac  delta  function  will  be  clear  from  the  context  in  which  it  is  used.  Our  main  result  con¬ 
cerns  the  perturbational  seismogram  Ss{c ,  f  ],  which  is  the  result  of  first-order  perturbation  of  the 
velocity  profile  and  source  wavelet: 

6s[c ,  f  j  =  8ut(z  =  0) 


where  the  perturbation  field  8u  satisfies 


P 


d28u 


V  •  8u 


8u 


o  2 

=  VV  •  8u  -  8p  ^4 
dt2 

=  ~Sf{t)8{x „  x2) 
=  0  t  «  0 


and  8p  =  —  2  8c  ■  c-3. 

Our  main  result  may  be  stated: 


The  velocity  profile  and  source  wavelet  perturbations  Sc,  8f  are  uniquely  determined  by  the 
perturbational  point-source  seismogram  8s,  under  hypotheses  Ml-MS  and  Sl-SS. 
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Note  that  we  do  not  assume  the  presence  of  an  isolated  reflector.  Thus  the  influences  of  Sf 
and  Sc  are  coexistent  in  time.  Nonetheless  they  can  be  recovered  separately.  In  our  view,  it  is 
this  feature  which  makes  the  present  result  interesting  despite  the  severe  restrictions  imposed  on 
our  model. 

This  result  concerns  the  perturbational  (“variable-background  Born  approximation”)  prob¬ 
lem.  The  technique  is  constructive,  and  estimates  Sc,  Sf  in  terms  of  St .  It  seems  clear  that  regu¬ 
larity  results  similar  to  those  explained  in  Symes  (1986b,  c)  could  be  combined  with  the  conclusions 
presented  here  to  yield  unique  determination  of  /  and  c  from  a .  Such  matters  will  be  discussed 
elsewhere. 

We  begin  the  study  of  this  problem  by  transforming  the  point-source  seismogram  to  a 
plane-wave  (p-tau)  section,  and  determining  its  structure,  in  Section  2.  We  use  this  structural 
information  in  Section  3  to  derive  a  constructive  procedure  for  separate  determination  of  source 
and  velocity  perturbations.  We  suggest  numerical  techniques  for  the  solution  of  least-squares  for¬ 
mulations  of  the  inverse  problems  and  present  the  results  of  numerical  experiments  based  on  one  of 
the  possible  implementations,  in  Section  4.  These  experiments  were  performed  by  Paul  Sacks  at 
Iowa  State  University.  Section  5  contains  the  proofs  of  the  main  mathematical  result  and  of 
several  technical  lemmas  needed  in  Section  3.  We  end  the  paper  with  a  discussion  of  possibilities 
for  relaxing  the  restrictions,  and  extending  the  scope,  of  our  results  (Section  6). 

2.  Structure  of  the  Perturbational  Seismogram. 

Our  main  tool  in  this  paper  is  the  plane-wave  decomposition:  since  the  coefficients  are 
independent  of  time  and  of  the  horizontal  coordinates,  the  Radon  integral 

//  dx  dy  u(x ,  y ,  z ,  t  +p  ■  x) 

produces  for  each  p  a  field  satisfying  a  system  of  partial  differential  equations  in  z,  t.  In  fact,  we 
are  only  interested  in  a  finite  depth  interval  0  <  z  <  zm„,  so  without  loss  of  generality  we  assume 


that 
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0<e.<e(z)<e*,  0  <  z 

for  suitable  constants  e*,e*.  In  Santosa  and  Symes  (1985,  Appendix)  it  is  shown  that  for 
sufficiently  small  p ,  the  plane-wave  component  of  normal  displacement 

U{z,r,p)  =  J  J  dx  dy  ut{x,y ,  z,r+  px) 


solves  the  boundary  value  problem 


(p-xp2) 


— — =  0 
dr2  dz 2 

-  x  4“  (°>  T>  P )  =  /  (r) 


(2.1) 


Urn  0  ,  r  «  0  . 

I _ J_ 

Recall  that  the  wave  velocity  is  given  by  e(z)  =  X2p  2(z).  After  normalization,  we  may  assume 
for  the  rest  of  the  paper  that  X  =  1,  so  that  the  mechanics  are  described  by  c  alone. 

The  equation  (2.1)  is  hyperbolic  only  so  long  as  c  |p  |  <  1  (precritical  slowness).  We  define 

zmtx(p)  =  wp{z-  c(z')p  <  1  for  0 <z'<z} 


and  consider  (2.1)  and  related  equations  only  in  the  slab  {(z,  t ):  0  <  z  <  zm«(P )}  for  each  p  - 

We  note  that  the  Radon  transform  integral  given  above  must  in  general  be  modified  by  the 
introduction  of  a  mute  or  cutoff  for  large  but  still  precritical  p  .  For  a  suitable  choice  of  mute,  the 
components  still  satisfy  the  plane-wave  equation  (2.1)  up  to  an  error  which  may  be  controlled  by 
the  techniques  presented  here:  see  Santosa  and  Symes  (1987). 

The  plane-wave  seismogram  (p-tau  section)  5  is  the  surface  (z  =  0)  time  history  of  the  plane 
wave  normal  displacement  fields 


S[c,f}{T,p)  =  U(0,T,p)  . 


The  topic  of  our  paper  is  the  “inverse”  problem: 


Given  a  measurement  of  G  of  plane-wave  reflection  data,  find  a  velocity- source  pair 
(c  ,  /  )  for  which 
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S[c,f}{r,p):  =  f/(0, 7-,p)  =  G(r,p)  . 

This  formulation  supposes  implicitly  that  the  (r,  p  )-pairs  at  which  the  data  G  is  given 
are  precritical  for  the  sought  velocity  profile  c ,  which  presumption  is  itself  a  nonlinear 
constraint  on  c . 

Now  it  is  natural  to  think  that  any  feasible  measurement  G  would  be  contaminated  by  error, 
hence  might  not  fit  any  model  exactly.  Accordingly,  it  is  popular  to  replace  the  “exact  inverse” 
problem  above  by  a  best-fit  formulation.  A  common  choice  of  error  measure  is  the  mean  square, 
which  leads  to  the  least  squares  inverse  problem: 

Find  e  ,  /  to  minimize 

JJ  dn(T,p)\S[c,f]{T,p)-  G(t,  p)  i2 

where  dft  is  a  measure  (i.e.  possibly  nonuniform  weight,  continuous  covariance  matrix). 

In  this  paper  we  shall  limit  our  weights  to  those  of  the  form 

</fl(r,p)  =  dr  dfi(p)  . 

Thus  time  points  are  uniformly  weighted  in  each  trace  (p  =  const.  ),  but  we  allow  the  traces  to  be 
weighted  according  to  the  measure  dp{p).  For  example,  if  dp(p)=dp,  then  all  traces  are  uni¬ 
formly  weighted  (over  the  domain  of  integration),  whereas  if  dp(p)  is  a  finite  sum  of  point  masses, 
then  the  error  measure  above  is  concentrated  on  the  corresponding  finite  set  of  traces. 

The  problem  stated  above  is  nonlinear.  To  begin  our  study  in  this  paper  we  shall  assume 
that  the  velocity  is  a  sum  of  a  (smooth,  slowly  varying)  background  velocity  c  and  a  (possibly 
rapidly  varying)  perturbation  Sc .  We  make  a  similar  assumption  concerning  the  source,  i.e.  that  it 
is  the  sum  of  a  background  wavelet  /  =5+  f  f  \  square-integrable,  and  a  perturbation  8f  .  We 
shall  study  the  corresponding  perturbation  SS  in  the  seismogram,  which  depends  linearly  on  Sc , 
6}  . 

Formally,  to  first  order,  the  result  of  perturbing  the  velocity  profile  c  *—  c  +  Sc  and  the 
source  wavelet  /  *—  f  +Sf  is  to  perturb  the  plane-wave  field  U  by  a  field  SU  satisfying 
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d25U 

dr2 

dSU 

dz 


d25U  25c  d2U 
dz2  c3  dr2 

{0,  r,  p )  =  <5/  (r) 


SU  =  0,  t  <  0  . 


(2.2) 


The  perturbational  seismogram  5S  is  the  surface  value  of  5U: 

SS(t,p)  =  6U( 0,  r,p). 

Presumably  S[c  +  5c,  f  +  5f  ]  « S[e  , / ]  +  5S .  Conditions  under  which  this  is  true,  i.e.  5S  is 
actually  the  derivative  of  S  are  discussed  in  Symes  (1986a). 

To  understand  the  structure  of  the  perturbational  seismogram  SS ,  we  now  introduce  asymp¬ 
totic  approximations  for  the  various  fields.  There  will  result  an  expansion  of  5S ,  which  is  of  cru¬ 
cial  importance  to  what  follows. 

As  mentioned  in  the  introduction,  we  consider  in  this  paper  only  deconvolvable  source 
wavelets  /  :  that  is,  we  assume  that  the  convolution  operator  with  kernel  /  has  a  bounded  inverse 
on  L2[ 0,  T]  for  any  T  >  0.  It  is  sufficient  (but  not  necessary)  that  the  Fourier  transform  /(w) 
have  a  uniformly  bounded  reciprocal.  (Of  course,  we  make  no  such  restrictions  on  5}  .)  Since  (with 
the  obvious  notation) 

S,  =  /  *  S5 

we  could  in  principle  apply  the  convolution  inverse  of  /  to  the  seismogram,  and  to  the  perturba¬ 
tional  seismogram.  This  amounts  to  replacing  5f  by  (f*)~l5f  ,  and  /  by  5,  in  all  of  the  above 
equations. 

The  field  U  is  now  the  plane-wave  impulse-response.  As  is  well-known,  this  field  possesses  a 
progressing  wave  expansion:  see  Courant  and  Hilbert  (1962,  Chapter  6)  for  generalities  and  Symes 
(1981,  Section  2)  for  computation  of  the  present  example.  We  obtain 

U{z,T,p)=  v{0,p)2v{z,p)2  H(t—  T{z,p))  +  R{z,T,p)  (2-3) 


where 
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'i  -c\z)r 


is  the  plane-wave  vertical  wave  velocity, 


is  the  vertical  travel  time  for  the  plane-wave  at  slowness  p ,  H  is  the  Heaviside  unit  step  function, 
and  R  is  a  continuous  remainder  term.  This  expansion  and  others  like  it  depend  for  their  validity 
on  the  existence  of  sufficiently  many  derivatives  of  e ,  which  we  assume.  See  the  concluding  section 
for  a  little  discussion  of  this  point. 

To  derive  the  promised  expansion  of  SS ,  first  consider  the  case  Sf  =0.  Then  a  Green’s  iden¬ 
tity  argument,  detailed  in  Santosa  and  Symes  (1987,  Section  5)  results  in  the  expansion,  expressed 
in  terms  of  7:  =  <51og  c  —8e  j c 

SS  |{/_o(r,p)  =  t»(0,  p){(l-  c2p2)~l7}  (Z(t,  p))  +  Kc~t  {t,  p)  (2.4) 

where  c0  =  c(0)  and  Ke  is  an  operator  of  Vol terra  type: 

Z[r,p) 

Ke7 (r,p)=  /  dzk'{p,T,z)  7(*) 


with  continuous  kernel  ke,  and  Z(r,p)  is  the  inverse  of  the  two-way  travel  time: 

*  =  Z{2T{z,p),p)  . 


On  the  other  hand,  a  perturbation  in  /  (=  6)  gives  the  solution  of  the  inhomogeneous  problem 


(_L_p2)iL_^L  w_0 

1  C2  P  ’  dT2  dz2 


=  -sf 


8U  a  0  ,  t  <  0 


which  is  in  turn  the  convolution  of  <5/  with  the  solution  of  the  same  problem,  but  with  Sf 
replaced  by  (Dirac)  6.  We  recognize  that  the  latter  problem  is  identical  to  (2.1)  with  /  =  —  S,  so 
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«$!&*)-«/  *u  |,_o  - 

Using  the  progressing  wave  expansion  we  see  that  this  is  proportional  to  the  integrated  wavelet 

f 

4>(t)  =  I  Sf 

0 

up  to  an  error  given  by  a  Volterra  operator  on  4>.  Introducing  the  velocity  trace 

V{T'P)  =  ■|^(0,rlp) 

we  can  re-write  the  above  relation  as 

(2-5) 

and  V(t,  p)=  v(0,  p)  S(t)  4-  (continuous  error)  is  an  invertible  (on  L2)  r-convolution  kernel,  known 
from  the  reference  seismogram. 

Combining  (2.5)  and  (2.4)  we  get 

SS(r,p)=  V*  4>(r,p)- v(0,p){(l-c2p2)-17}(Z(r,p))  +  K7(r,p)  .  (2.6) 

Observe  that  V*  has  a  simple  convolution  inverse 

vi{r,P )  =  v(0,p)_l5(r)  +  *1(T,p), 

where  is  continuous  and  causal  for  each  p  .  We  shall  use  this  fact  in  the  algorithm  in  Section  4. 

3.  Separation  of  Wavelet  and  Velocity  Perturbations 

It  is  evident  from  (2.6)  that  the  p -dependence  (“moveout”  and  scaling)  of  Sf  and  Sc ,  as  they 
appear  in  SS ,  are  different.  In  this  section,  we  draw  a  consequence:  that  both  Sf  and  Sc  may  be 
determined  from  SS . 

We  do  not  address  the  quality  of  this  determination  here,  or  even  the  precise  conditions  under 
which  it  holds.  We  want  the  reader  to  be  aware  of  the  limitations  of  such  formal  analysis:  it 
gives  no  insight  whatsoever  into  the  effectiveness  of  any  algorithms  developed  to  exploit  the 
theoretical  observation.  Such  insight  can  only  be  gained  by  a  rigorous  (i.e.  correct  and  complete) 
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analysis,  which  we  give  in  Section  5. 

In  keeping  with  the  formal  approach  of  this  section,  we  will  drop  all  “Volterra”  terms  (i.e.  K 
in  (2.6))  as  negligible.  The  sense  in  which  this  is  actually  true  is  explained  in  Section  5. 

Thus  (2.6)  becomes 

SS  »  «(0 ,  p)  (^(f)-  [(1-  cvr^K^np)))  • 

(Here  we  have  used  the  remark  at  the  end  of  Section  2,  writing  V(r,  p )  ~  v(0,  p )  6(z). 

Thus,  except  for  scale,  4>  appears  without  any  p -dependence  at  all.  Accordingly,  form  the 
auxiliary  data  set 

G{t,  Pi,p2)  =  «(0,  Pi)~lSS  {T,pl)-v(Q,p2)~l6S(T,p2) 

=  [(1-c2P2  )~1i](Z{t,p2))-  [(l-c2p12)_17j(Z(r,  Pl))  • 

Suppose  that  pt  >  p2.  Then  v(z  -Pi)  >  v(z,  p2)  for  all  z,  so  the  prwave  reaches  any  given  depth 
earlier  than  the  prwave.  Accordingly, 

°(*.Pi>P2):  =  Z(2T(ztpi),p2)<z  .  (3.1) 

The  transformed  auxiliary  data  set 

G{z,p i,  p2):  =  (1  -  c2{z)pf)G(2T(z,pi),pu  p2) 

is  the  left-hand  side  in  the  functional  equation  for  7: 

G(z ,Pi,  p2)  =  l(z)  +  P(z ,  p x,  p2)  l{ot(z ,  Pi,  P2))  (3-2) 

where 

^(z,Pi,p2)  =  (1-  c2(r)Pi)(l-  c2(Q(«,Pi,P2))2P22r1  • 

Because  of  the  “delay”  inequality  (3.1),  the  equation  (3.2)  gives  7(2),  for  any  z,  as  a  combination 
of  the  data  G  and  the  value  of  7  at  some  shallower  z'(<  z).  Thus  it  would  seem  reasonable  that  7 
would  be  entirely  determined  if  we  know  7(z)  for  some  very  shallow  near  surface  layer  0<z  <  z., 
for  then  we  could  work  our  way  downward,  using  (3.2)  recursively,  to  obtain  7  at  any  depth. 
Thus  some  restriction  on  the  behaviour  of  7  near  z  =  0  would  allow  us  to  determine  7. 
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While  this  is  true,  it  is  unnecessary  to  restrict  7  near  the  surface,  because  of  the  properties  of 
the  coefficients  a  and  0.  This  result  is  somewhat  technical,  and  will  be  given  in  detail  in  Section  5. 

In  sum,  the  functional  equation  (3.2)  determines  7(2)  (without  any  restrictions).  Then  7(2) 
may  be  substituted  in  (2.6),  which  may  then  be  solved  for  ^(r)  (for  any  p).  Thus  both  7(2)  and 
0(r)  are  determined. 

In  the  next  section,  we  consider  the  algorithmic  implications  of  this  observation. 


4.  Numerical  Solution  of  the  Inverse  Problem. 

In  this  section  we  show  how  the  considerations  of  Section  3  yield  an  iterative  algorithm  for 
solution  of  the  least-squares  inverse  problem  in  the  (r,p)  domain,  and  discuss  its  implementation 
for  a  special  case. 

Thus  we  attempt  to  minimize,  for  a  given  (p ,  r)  data  set  g , 

^nrnxfl’) 

\\S(c,f)-g  \\*=  J  dfi(p)  J  dr\S-g\2  (4.1) 

0 

in  the  notation  of  Section  3.  In  fact,  we  assume  here  that  on  Tmax(p )  =  I  ssd  that  the  domain  of 
integration  consists  of  precritical  (r,  p )  pairs  for  all  velocity  profiles  to  be  constructed  during  the 
iteration.  This  assumption  a  priori  restricts  the  velocity  profiles  and  allows  us  to  use  only  a 
minimal  set  of  “safely  precritical”  data.  Obviously  relaxation  of  these  restrictions  would  be  desir¬ 
able.  For  some  discussion  see  Santosa  and  Symes  (1987),  Sections  2  and  7. 

From  (2.6)  we  can  write 

DS  { e  ,  $](5c ,  6f  )  =  E  {<f>,  7)  +  K  (0, 7)  (4.2) 

where 


DS  [c,S}{8f  ,6c)  =  6S 


is  the  directional  derivative  of  the  seismogram  5, 


E(*,7)(r,p)  =  v  (0,  p) 


•  > 
<?Kt)+  [(1  —  c2P2)-1-7  {Z{p,t)) 
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and  K  is  a  Volterra-type  integral  operator  with  continuous  kernel. 

We  shall  suppose  that  the  wavelet  is  a  slowly-varying  perturbation  of  6: 

f(r)  =  S(z)  +  f  i(z) 

where  f  i(z)  is  square-integrable,  as  mentioned  in  the  introduction.  Thus  the  integrated  wavelet  <f> 
is  continuous  for  r>  0,  and  satisfies 

lim  <I>(t)  =  1  . 

r— 0+ 

Then,  since 

DS[c,f  ](&  ,Sf)  =  f*DS[c,  6}(6c  ) 

in  fact  DS[c ,  f  ](6c  ,Sf )  is  given  by  the  same  expression  (4.2),  where  K  is  modified  to  include  con¬ 
volution  by  the  square-integrable  part  of  /  . 

The  Gauss-Newton  method  for  minimizing  J(c,f)  updates  a  current  estimate  (ce,fe)  by 
solving  a  linear  least-squares  problem: 

e+—  cc  +  6c  ,  /+  =  fe  +Sf 

where  (Sc  ,  <5/  )  minimizes 

I  \DS[et ,  fc  }(Sc  ,  6f)-(g-S[ce,  fe  ])  ||»  .  (4.3) 

The  minimizer  of  (4.3)  satisfies  the  normal  equations 

DS  [ce ,  fe]'DS  [ec ,  { c\(Sc  ,  Sf  ) 

—  DS  (cei  fe]'(g  —  S  (ce,  f  })  (4'4) 

where  DS '  is  the  adjoint  of  DS  (regarded  as  a  linear  map  acting  on  (4>,  q))  from 
L2[0,  T]  X  L2[0,  Z]  —  L2([0,  T j  X  supp  d/t,  d  rd/t(p)). 

Since  K  in  (4.2)  is  Volterra,  it  represents  a  small  perturbation  if  Z,  T  are  small  enough. 
This  suggests  replacing  DS  and  DS *  in  (4.4)  with  the  simpler  operators  E  and  E*,  with  the  latter 
given  by 
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e'$  =  (e;$,  e7*$) 

E/$(r)  =  J  dfi(p)  v  (0,p ) $(r,  p ) 

e;^)  =  2  fjM(p)  Mo.p)y2^,p),p)  , 

«(z,p)(l-c2(i)p2) 


It  follows  from  the  heuristics  of  Section  3  and  the  estimates  in  Section  5  that  E*E  is  inverti¬ 
ble  —  in  fact  positive  definite  —  and  E  and  E*  are  easy  to  evaluate.  Thus  the  modified  normal 
equations 

E*E(<6,q)  =  E  ‘{g-S[ce,fc\)  (4.5) 

are  uniquely  solvable,  and  an  iterative  method  like  conjugate  gradients  should  be  quite  efficient. 

In  rough  outline,  the  resulting  algorithm  is:  repeat  until  convergence 

(1)  compute  the  seismogram  S  [c  ,  /  ]; 

(2)  solve  the  modified  normal  equations  (4.5)  for  <t>,  7; 

(3)  update  c  ♦—  ee7,  /  «—  /  +  <£';  go  to  (1) 


Essentially  the  same  approach  (replace  DS  by  E)  is  used  with  success  by  Sacks  and  Santosa 
(1987)  in  recovering  c  alone  (they  consider  the  “consistent”  case  and  solve  a  functional  equation 
but  the  idea  is  very  similar). 

Since  the  replacement  DS  ■*—  E  involves  only  a  small  change  for  r  and  z  small,  we  would 
expect  this  algorithm  to  behave  like  the  Gauss-Newton  iteration  near  r  =  z  =  0.  The  estimates  of 
Section  3,  which  show  that  DS  is  bounded  below,  would  suffice  to  show  that  the  Gauss-Newton 
algorithm  is  convergent  if  5  were  differentiable.  Unfortunately,  differentiability  of  S  requires  more 
constraints  on  c  than  we  have  included  here  —  essentially,  c  must  be  relatively  smooth.  The  rea¬ 
sons  for  this,  and  a  convergence-inducing  regularization  of  the  least-squares  problem,  are  discussed 
in  Symes  (1986c).  Granted  that  this  regularization  has  been  implemented,  either  by  restricting  the 
class  of  admissible  c ’s  or  by  adding  a  penalty  term  to  J,  Gauss-Newton  will  converge  and  thus  so 
will  the  modified  algorithm  given  above. 
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We  will  now  consider  a  special  case  of  this  procedure.  More  general  examples  will  be 
explored  elsewhere. 

We  take  (for  Pi  >  p2>0) 

dfi(p)  =  S(p  —  Pi)  +  S(p  - p2 ) 

which  amounts  to  choosing  for  our  data  exactly  two  plane-wave  traces.  Now  we  should  have  “just 
enough”  data  to  solve  the  problem:  two  time  series  to  determine  two  time  series.  In  fact,  a  minor 
modification  of  the  discussion  of  Section  3  leads  to  the  conclusion  that,  for  appropriate  z  and  r 
intervals,  the  mapping 

<t>  E(*,7)  (p1(  ) 

7  J1  “  (e(<M)  (?2,  ■). 

is  invertible  (as  is  the  analogous  map  with  E  replaced  by  DS[e ,  /  ]).  It  follows  that  E*  is  inverti¬ 
ble  as  well,  so  the  modified  normal  equations  (4.5)  are  (in  this  special  case)  equivalent  to 

E(^,7)  =  9  -S  [e,f]  (4.6) 

which  is  an  approximation  to 

DS[e,f}(Se,6f)=  g-S[e,f]  .  (4.7a) 

Now  (4.7)  defines  the  Newton  update  for  the  functional  equation 

S[c,f)  =  g  (4.7b) 

so  we  would  expect  iterative  use  of  (4.6)  to  yield  an  approximate  solution  of  this  functional  equa¬ 
tion,  rather  than  merely  a  least-squares  solution.  This  is  precisely  the  analogue  of  the  algorithm 
presented  in  Sacks  and  Santosa  (1986). 

The  prescription  for  solving  (4.6)  is  actually  given  in  Sections  3  and  5;  in  outline,  repeat  until 
convergence  the  following  steps: 

(i)  Create  the  auxiliary  data  set 


G(r)  =  Vx*  {g  -  S  [e  ,  f  ])(r,  Pl)  -Vl*(g-S[c,f  ])(r,  p2) 


(4.8a) 


IS 


where  Vj  is  the  convolution  inverse  of  the  (current)  velocity  trace; 

(ii)  Create  the  auxiliary  data  set 

G(z)  =  (l-c(zfp?)G(2T(z,Pl)); 

(iii)  Solve  for  7  the  functional  equation 

G(z)  =  7 (*)  +  £(*)  7  («(*))  (I+H^7(jr)  (4.8b) 

where 


a(z)  =  Z(2T(z,Pl),p2)  (4.9) 

1-  c(z)2Pi 

l-c2(a(z))pi 

using  the  Neumann  series 

OO 

(I+TT)-1  EC-TT)*  (4.10) 

the  convergence  of  which  is  guaranteed  by  the  properties  of  a  and  /?  and  by  Lemma  1 
of  Section  5. 

(iv)  Compute  <t>  from  the  relation 

V\  *  j(?  ~S  [c,/])(r,  pi)  v  (0,  p)  [(1  -  c2p12)_11f](^ (r,Pi)) 

=  Hr)  ; 

(v)  Update  c  and  /  : 


c 


ce 


1 


f  «-/+*'• 


Remark.  These  steps  may  be  simplified  further  by  noticing  that 
U1(r,p)  =  w(0,  p)-1'5(r)  +  jfcj(r,  p),  as  noted  earlier,  so  that  the  convolution  V,  *  may  be  replaced 
by  multiplication  with  v(0,  p)-1  at  the  cost  of  another  (“Volterra”)  error  of  the  type  we  are 
already  ignoring.  This  was  also  done  in  Section  3. 

We  have  implemented  a  single  step  of  this  procedure.  We  compute  the  plane-wave  seismo¬ 
gram  S  [e,f\  using  a  standard  leapfrog  finite-difference  method,  applied  to  the  first-order  system 
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equivalent  to  the  wave  equation,  in  which  the  velocity  is  one  of  the  dependent  variables.  These  and 
other  functions  were  represented  by  grid-functions  on  suitable  space-time  grids. 

To  obtain  the  unknown  <5c(z)  we  must  compute  the  function  G(z),  defined  in  (4.8a,  b),  then 
evaluate  the  series  (4.10).  The  compositions  with  travel-time,  and  with  the  depth  change-of- 
coordinates  (4.9),  are  computed  using  piecewise  linear  interpolation,  the  accuracy  of  which  is  com¬ 
patible  with  that  of  the  difference  scheme.  We  replaced  convolution  with  V  by  its  discrete  version: 

(V  *4>)(r,)~  £  At 

* 

and  deconvolution  by  V  (i.e.  convolution  with  V\)  with  the  solution  of  this  triangular  system  by 
back-substitution. 

A  synthetic  example  of  this  procedure  is  presented  in  Figures  1-5.  Figure  1  shows  the  refer¬ 
ence  velocity  c  and  the  velocity  to  be  recovered  c  +  <5e ,  and  Figure  2  shows  the  corresponding 
reference  and  perturbed  wavelets,  for  the  two  slownesses  we  chose  P2  =  0  and  p!  =  .4.  In  Figure  3 
we  show  S  [e,  /](r,  p,- ),  *  =  1,2  and  g{r,  p,  )  =  S  [c  +  Sc ,  f  +<5/](r,  pt  ),  «'  =  1,2. 

The  velocity  perturbation  Sc  is  computed  as  discussed  above;  it  was  necessary  to  use  three 
terms  in  the  series  (4.10).  The  resulting  estimate  of  c  +  Sc  is  displayed  in  Figure  4,  along  with  the 
exact  velocity  profile.  Finally,  the  exact  and  recovered  wavelets  are  shown  in  Figure  5. 

We  emphasize  that  this  example  illustrates  our  contention  that  the  source  wavelet  (perturba¬ 
tion)  may  be  recovered  even  without  the  presence  of  a  strong,  clean  reflection  in  an  otherwise  quiet 
part  of  the  section.  In  this  example,  the  perturbations  Sc  and  Sf  yield  trace  perturbations  Ss 
which  are  completely  time-coincident. 

Note  that  in  the  interesting  special  case  that  c(z)  =  1,  f(t)  =  S(t),  the  above  procedure 
simplifies  substantially,  since  5  [c  ,  /  ]  can  then  be  computed  exactly. 

Since  this  procedure  approximates  one  step  of  Newton’s  method  for  the  functional  equation 
(4.7b),  we  would  expect  the  accuracy  to  degrade  as  Sc  becomes  larger.  To  illustrate  this  degrada¬ 
tion,  we  show  in  Figure  6  the  perturbed  surface  traces  S  [c  +  Sc ,  f  +  Sf  ]  for  p2  =  0,  pt  =  .4  and  for 
Sc  and  Sf  representing  1%  of  the  “energy”  (L2-norm)  of  c  and  / .  The  recovered  wavelets  and 
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velocity  are  compared  with  the  targets  in  Figures  7  and  8.  Figures  9,  10,  and  11  repeat  this  com¬ 
parison  for  10%  perturbations  6c  and  6f  . 

5.  Estimation  of  the  Wavelet  and  Velocity  Perturbations. 

Our  aim  in  this  section  is  to  show  that  the  differential  section  <55  dominates  the  perturba¬ 
tions  q  =  61og  c  and  <j>  =  j  6f  : 

Ik  M2+  IMi2  <  C  ||55  ||2 

where  the  vertical  bars  denote  appropriate  L2-norms.  Such  estimates  are  necessary  to  ensure 
stable  linear  estimation  of  6c  and  6f ,  and  also  form  the  heart  of  the  construction  of  nonlinear 
least-squares  estimators  for  e  and  /  . 

Our  derivation  will  be  constructive,  and  will  justify  a  computational  technique  which  we 
explored  in  the  last  section,  as  well  as  the  heuristics  of  Section  3. 

Recall  that 

is  the  convolutional  inverse  of  V(r,  p),  with  k\  continuous  and  causal.  From  the  structure  of  Vx, 
and  the  expansion  for  <55  in  (2.6),  the  auxiliary  data  (Section  3)  must  have  the  form 

G{r,  Pi,  Pz)  =  (l-c2Pi)_17(^(r,p1)) 

-  (1  -  c2P2)_!7  {Z{t,  p2))  +  {t,  V  i,  P2) , 

where  is  an  integral  operator  with  piecewise  continuous  kernel  k i(pi,  p^,  r,  z)  supported  in 
(0  <  z  <  max  ( Z(r ,  p !),  Z{r,  p2))}.  This  last  term  was  dropped  in  Section  3. 

We  have  already  indicated  that  when  p{  >p2, 

a  [z,  Pi,  P2):  =  %  (2T(r,  pi),  p2)  <  r  . 


We  will  now  quantify  this  inequality. 
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The  quantity  a  (z ,  p  ,,  p 2)  is  crucial  to  the  argument  which  follows:  it  gives  the  depth 
reached  by  a  wave  traveling  at  the  slower  wave  speed  v(- ,  p2)  in  the  two-way  time  taken  by  a 
wave  traveling  at  the  faster  wave  speed  v( • ,  pt)  to  reach  depth  z.  thus  it  represents  the  “spatial 
delay”  caused  by  replacing  the  (faster)  t>(- ,  Pl)  by  the  (slower)  v( ■ ,  p2).  We  shall  assume  that  we 
consider  only  pairs  (r,  p )  for  which  2  =  Z(r,  p )  satisfies 


-  p2  >  A2  for  0  <  z'  <  z 

for  some  fixed  A  >  0.  Denote  by  zm„(p)  the  largest  z  for  which  (5.2)  holds. 

To  take  into  account  the  possibly  finite  duration  of  the  p-r  traces  (0  <  r  <  Tm%x(p )),  set 


(5.2) 


zmtL x[P  )  —  (z  majt(P  )»  %(  T msj((p  ),  p  ))  • 

(Note  that  z  m„(p )  =  00,  for  instance,  for  p  =  0,  if  e  is  bounded  on  0  <  z  <00  and  Tm„(0)  =  00). 
Correspondingly,  set 


rmM  (P)-2T(  (p)>p)  • 

Then  for  0  <  z  <  zm„(p )  a  little  algebra  yields 

t(z,p2)>(i  +  ^-(p?-p2))T(z,Pi), 

whence 


z  -  Z{2T  (z ,  p  1),  p2)  >  c.  T  (z ,  p2)  —  T  (z,pi) 


>  -P2  )  TiZ’Pl) 

>  —T~(P\  ~Pi)z  ■ 


Thus 


a  {z  t  P  U  P  2)  <  {Pl,P2)z 


(5.3) 


with  a  '  =  1  — —  c.3A(p2  —  p2  )  <  1.  Similarly,  it  is  easy  to  see  that  a'>a,  >0  for  some  a. 
4 


independent  of  z ,  p . 


20 


Recall  the  definition  of  G(z,  pb  p2)  and  its  representation  (3.2): 

G  (*,Pi,pa):  =  (1  -  c2(z)p?)G{2T  (z,  px),  pu  p2)  ■ 

(5.4) 

=  7 (z)+0(z,pl,p2)'y(a(z,Pi,p2))  +  K2i(z,p1,p2 ) 

where  K2  is  another  integral  operator  with  piecewise  continuous  kernel  k^p  p2,  z ,  z')  supported  in 
(0 <  z'<  z},  and 


Hz>PuP2)  =  (k  -  c2(z)pl){l  -  e2{a{z ,  pu  p2))  pi)  1  . 


(Thus  (5.4)  is  the  precise  statement  of  which  (3.2)  is  an  approximation.)  Note  that 


0(0-  Pl,P2)  = 


1-  c2(0)pf 

1  -  c2(0)p| 


<  1  . 


(5.5) 


Now  it  follows  from  (5.3),  (5.4),  (5.5)  and  Lemma  3  below  that 

llT,NL2[0,ImJ(,1)|  ^  C'  NG(^Pl.P2)lli2|0tWp^ 

<C2  ||G(  ■  ,  P,,  P2)  I|i2r0r  If  )| 

(5.6) 

<^{ll^(-,P1)lli2,o,WPl)1 

+  ll55(--P2)lli2,0,Wp2)]}  • 

From  the  foregoing  calculations  and  the  statement  of  Lemma  3,  the  constants  Cu  C2,  C3  evi¬ 
dently  depend  on  e.,  A,  and  p  2  —  p2  >0. 

To  express  the  dependence  of  6c  and  6f  on  the  entire  precritical  p-tau  section,  choose  a  posi¬ 
tive  measure  dp.  and  set 


11-55  ||^  =  /  dp(p)  J  dr  |55(r,p)  |2  . 
o  o 


Thus  ||  || 2  is  a  weighted  mean-square  error  measure.  For  example,  if  we  choose  dp(p)=dp 

(Lebesque  measure),  the  all  p-tau  traces  are  weighted  uniformly,  whereas  choosing  dp(p)  to  be  a 
linear  combination  of  point  masses  has  the  effect  of  selecting  a  discrete  set  of  slownesses. 

Now  suppose  that  for  some  e  >  0, 

{(pi,  p2):  I P i  —  P2  I  >f}D  (s«PP  dpXsvpp  dp)  5*  0  . 
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Let  pmin  =  inf{  ]p  |:  p ,  p  4-  £  £  aupp  /*}.  Then  from  (3.6)  follows  immediately  the  estimate 

m“lrm«(Pl)'rintt(P2)l 

<Gi  /  dp(Pl)  J  d^{p2)  Jdr  \G{t,PuP2)\2 

Ipi-p2I>«  o  l15'! 

<^6  \\6S  ||* 

where  C4  and  Cg  depend  on  G  and  n  as  well. 

Perhaps  a  refined  analysis  of  these  constants  is  possible,  in  the  spirit  of  Santosa  and  Symes 
(1987). 

Finally,  with  7  known,  <t>  can  be  extracted  from  any  of  the  components  via  (2.6),  whence  esti¬ 
mates  for  4>  follow  immediately. 

To  complete  this  section  we  give  estimates  needed  above,  which  establish  the  invertibility  on 
L2  of  operators  of  the  form 

I  +TT+  K 

where TTis  multiplicative-delay:  - 

TT«(i)  =  P{z)u{oc{z)), 

with  a(z)  <a*<l,0<z<l  and  K  is  an  integral  operator  of  Volterra  type 

2 

Ku(z)  =  /  . 

0 

For  convenience  we  restrict  our  attention  to  the  scalar  case.  Analogous  results  for  vector-valued 
functions  follow  from  trivial  modifications  of  the  proofs.  We  begin  with  the  case  k  =0. 

Lemma  1 .  Suppose  that  a,  PE  C^O,  l]  satisfy 


(1)  or(0)  =  0 

(2)  for  some  a.,  a’  with  0<a.,a*  and  a  *  <  1 , 
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l 

(3)  /9(0)(a'(0))  '  <  1. 


a'(z)  >  a» 
a(z)  <0**2, 


0  <  2  <  1 


Then  for  any  fi  <  1,  the  operatorTE  Lp{ 0,  l]  — *•  Lv[ 0,  l]  defined  by 


Th(2)  =  P{z)  “(<*(2)) 


satisfies 


IP"  II,  <#* 

for  some  N  depending  on  ||a||c,(01),  \\p  ||cl|01j,  a . ,  a  * ,  £(0)(a'(0))  p ,  and /i. 


Proof.  Define  sequences  an,Pn  in  C71  [0,  l]  by 


Then 


ao(2 )  =  2 

<*1(2)  =  0(2) 

“.(2)  =  a(a._i(2)) 


for  0  <  z  <  1, 


n  >  1 


«_,(?)  =  a  ‘(z)  for  0  <  2  <  a(l) 

a-n{z)  —  Q-n+1(Q-l('Z))  for  0  <  2  <  Of»(l)  ,  «>1 


0l(*)  =  P(z) 

P*(z)  =  0(*)0»-l(«(*)) 


for  0  <  z  <  1 ,  n>l 


=  II  0(or*(*)) 


TT«(z)  =  /?,(*)  «(“■*(*))  . 


Note  that 


0,(2)  <  -1(2)  <  •  •  <(q*)*2,  0<z  <  1 

so  that  aB(l)  <  (a  *)n  — *■  0  as  n  —*•00.  Choose  5  so  that  for  0  <  z  <  S, 

P(z)(a'(z))~  <  l[l  +  /?(0)(a'(0)f ']:  =  X  <  1 


and  select  n  so  that  (a*)"  <5. 


23 


Note  that 


<*n  =  <*  '  a»-i  •  a,_i 

*  —  1 

=  ....=  FI  a'  •  a* 

km C 


Note  also  that  if  z  £  [0,  at(l)j,  then 


Thus 


Set 


ay  •  a_*(z)  =  ay_t(z) ,  >  =0,1,2, 


1 

ll^tt  lip  =  S  d*  \Pn*  ’aN  |p 
0 

^/vO) 


*w(>) 

=  / 


r  I/»n  •«-*!'  ,  „ 
=  J  ——r — —  l«  r 

o  aN  a_N 

l«  lP 


\fi  ‘  Ctk-N  \P 
kmc  a'  •  ak_N 


(5.8) 


Then 


<?:-  ||/J{oO  Ml 


"-1  \P  •  ak-N  V 

k- o  a  *  ak_N 


<  C> 


(5.9) 


on  [0,  aN(l)].  Also  note  that,  if  z  G  [0,  aN(l)]  and  k  >n  ,  at_w(z) =  ak^N(aN(w))  =  ctk(w)<  5  for 
some  w  £  [0,  l].  Thus 

\P  *«*-N  I” 


„„  •(*)<*' 

“  '  <Xk-N 

for  k  >  n,  0  <  z  <=  aiV(l)  . 


(5.10) 


Taken  together,  (5.8),  (5.9),  and  (5.10)  imply  that 


IP?'  ||;  <  Cr\^-’t)p  1 1 u  Wg 


for  N  >  n .  If 
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»r  -  ..  .  log  C  —  log  ft 
iv  >  n  +  — — - — — — 

(—  log  X) 

then 

IP*'.  iu  <  ^  lit.  II, 

as  required.  Clearly  C ,  \  and  n  depend  on  the  stated  quantities. 

We  now  consider  operators  of  “Volterra”  form,  that  is, 

l 

K«(z)  =  /  dgk(z,g)u{g)  . 
o 

As  is  well  known,  for  continuous  kernels  k,  such  operators  have  spectral  radius  zero  as  operators 
on  Lp .  See  e.g.  Widom  (1969),  p.  9-11. 

Actually,  slightly  less  than  continuity  is  required.  We  consider  operators  K  defined  by  ker¬ 
nels  k ,  possibly  matrix-valued,  for  which  the  function 

*(*)-*(•,*) 

lies  in  C°([0, 1],  L2[0, 1])  with  the  property  that  supp  Ic  (z)C  [z,  1],  0<  z  <  1. 

Note  that  the  Volterra  property  is  equivalent  to  the  statement  that 
supp  «  c  [*i,  l]  =>  svpp  Kt»  C  [*i,  1]  ■ 

Lemma  2.  Suppose  that  k  G  C°( [0,  l] ;  L2[ 0,  l])  satisfies  supp  Ic  (z)  C  [z,  l]  and  define  for  u  G  L2[ 0,  l] 

l 

K«  =  J  dz  f  (z)  u(z)  G  L2[ 0, 1]  . 
o 

Then  for  any  positive  /i  <  1  there  exists  a  positive  integer  N  depending  on 


and  on  n,  for  which 

I|K"||2<1. 


Proof.  (This  is  a  small  modification  of  standard  argument.  See  e.g.  Widom  (1969).)  Define 


the  iterated  kernels  k  n  by 
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k  i  =  Ic 

=  J  dz  k'H^(z)k(zh  z)  n  =  2,3,4  •  •  ■ 

One  easily  verifies  that  £  ^([O,  l];L2[0, 1]).  Because  of  the  observation  preceding  the  state¬ 
ment  of  the  lemma,  supp  Ic  »(z)C[z,  l]  as  well. 

Note  that 

ll^2(2i)  II  =  J  dz  Ic  i(z)*(*>*i)  || 

1 

=  ||/<fc*~(2)*(z,Zi)|| 

*1 

1 

<  J  dz  ||*'(z)  ||  |*(z,zi)  | 

1 

<Cjdz  |*(z,z,)| 

<C\  l-z,)2 


where  we  have  written  G=  sup  ||Jt'(z)||.  Assume  that  ||fy(z)  II  <  — --(1-  ^ — , 

(n\f 

j  =1,2,...,  n.  In  general, 


ll^»+i(*i)  II  <  J  dz  ||f.(OII  l*(*>zi)l 

*1 

/-m+l  1  JL 

<  - -  f  dz  (1- z)2  \k{z,Zy)\ 


< 


(n!)2  ‘l 
(7*+1 


fdz(l-z)' 


(n!)2  1'* 

0n+2 


ll*(*l)ll 


rU-*i) 


n+j 


((n  +  l)!)5 


Thus  ||KS  ||  <  1  for  large  n,  as  desired. 


q.e.d. 


Lemma  3 
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Suppose  that  k  is  a  Volterra  kernel  as  described  in  Lemma  2,  and TT  i«  a  multiplicative  delay 
operator  defined  by  scalars  a,  /3  obeying  the  conditions  of  Lemma  1.  then 

I  +TT+  K 

is  invertible:  in  particular,  for  some  £  >  0  depending  on  the  quantities  described  in  Lemmas  1  and 

5, 

11(1  +T+  K)«  ||  >e  ||«  || 

for  u  C£2[0, 1]. 

Proof.  The  convergence  of  the  Neumann  series  for  I-fTTin  operator  norm  is  a  consequence  of 
Lemma  1.  Thus  (H+TI)_1K  is  a  bounded  operator  on  L2[0,  lj,  given  by  a  kernel 

fi(*)  =  (i+uTMO  • 

The  Volterra  property  of  k j  follows  from  the  delay  property  (a  <  1)  ofH  and  an  estimate  for  kl 
in  C'°([0, lj;  L2[0,  lj)  follows  directly  from  Lemma  1  as  well.  Thus 

||(i+Tr+K)«  II  =  ||(i+'H)(n  +  (i+Ti)-1K)«)|| 

>£1||(I  +  (I+'Il)-,K)ti  || 

><  ||«  II 

by  the  application  of  Lemma  2.  q.e.d. 

6.  Discussion. 

In  this  section  we  discuss  briefly  the  consequences  of  weakening  the  hypotheses  M1-M3,  S1-S3 
stated  in  Section  1,  which  underlay  our  arguments  but  which  are  unacceptably  restrictive  for  prac¬ 
tical  application. 

The  extension  of  the  perturbational  results  to  nonlayered  background  media  is  straightfor¬ 
ward,  to  some  extent,  so  long  as  the  background  is  smooth.  In  some  sense  the  literature  on  migra¬ 
tion  is  concerned  exactly  with  this  extension:  for  some  recent,  mathematically  correct  results  for 
the  fixed  /  case  see  Beylkin  (1985)  and  Rakesh  (1986).  Results  analogous  to  those  presented  here 
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should  be  achievable  in  this  context.  The  (full)  nonlinear  nonlayered  problem  is  much  more 
difficult;  for  some  limited  results  see  Sacks  and  Symes  (1985),  Symes  (1986a). 

The  isotropic  elasticity  model  is  the  general  choice  for  the  “premium”  model  level;  see  e.g. 
Tarantola  (1984).  Some  results  on  the  determination  of  a  layered  elastic  medium  from  various 
point-source  data  sets  may  be  found  in  Clarke  (1984),  Yagle  and  Levy  (1986),  Carazzone  (1986), 
and  Sacks  and  Symes  (1987).  On  the  other  hand,  at  least  transverse,  and  possibly  more  general, 
anisotropy  is  evident  in  reflection  data  in  some  locales  (Thomsen  (1986));  Coen  and  Meadows 
(1985)  give  some  results  concerning  anisotropic  elastic  inversion.  We  see  no  difficulty  in  extending 
the  results  reported  here  to  these  more  general  layered  models. 

The  smoothness  question  is  quite  interesting  and  is  presently  inadequately  understood. 
Symes  (1986c)  gives  a  review  of  the  extent  to  which  smoothness  constraints  can  be  relaxed  for  a 
single-plane-wave,  layered  medium  problem,  and  Sacks  and  Symes  (1985)  gives  some  idea  of  the 
difficulties  which  arise  in  the  study  of  several-dimensional  problems  in  the  presence  of  limited 
smoothness.  An  analytic  approach  to  the  nonlinear  problem  requires  that  this  issue  be  addressed, 
as  perturbations  and  backgrounds  must  then  be  regarded  as  having  the  same  degree  of  smoothness. 
On  the  other  hand,  nonsmoothness  may  have  very  favorable  consequences^  as  will  be  discussed 
below  in  relation  to  band-extrapolation,  and  must  certainly  be  regarded  as  a  feature  of  the  actual 
parameter  distributions  in  earth  materials. 

As  mentioned  in  the  introduction,  the  point-source  assumption  seems  adequate  for  reflection 
seismology  on  physical  grounds.  On  the  other  hand,  common  land  and  marine  energy  sources  are 
strongly  anisotropic.  To  some  extent  this  feature  might  be  modeled  by  a  multipole  source  term, 
and  the  moment  tensor  components  recovered  by  a  generalization  of  our  technique.  This  seems  an 
important  matter  for  further  work. 

Perhaps  the  most  unrealistic  assumption  of  our  work  is  S3  (quasi-impulsive  nature  of 
sources):  it  is  simply  violently  wrong.  Typical  reflection  seismic  sources  have  significant  energy 
content  in  the  range  4-60  Hz,  at  best.  Thus  a  reasonable  isotropic  point-source  model  should 
involve  a  very  non-impulsive  wavelet. 
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In  the  last  several  years,  considerable  numerical  evidence  has  emerged  to  support  the  conten¬ 
tion  that  bandlimited  point-source  data  does  determine  layered  velocity  profiles  by  least-squares 
data  fitting:  see  for  example  McAulay  (1985)  and  Kolb,  Collino,  and  Lailly  (1986).  Santosa  and 
Symes  (1986)  give  a  comprehensive  analysis  of  this  problem,  and  show  that  success  of  the  least- 
squares  approach  is  crucially  dependent  on 

(i)  sufficient  (precritical)  aperture  and 

(ii)  sufficient  reflector  density,  i.e.  sufficient  lack  of  smoothness. 

Thus  the  nonsmooth  distribution  of  parameters  may  have  as  a  positive  consequence  the  feasibility 
of  seismic  inversion. 

We  believe  that  the  ideas  of  this  paper  and  those  of  Santosa  and  Symes  (1986)  might  be  syn¬ 
thesized  to  obtain  codetermination  of  velocity  profiles  and  bandlimited  point-sources,  under 
appropriate  conditions.  In  support  of  this  conjecture  we  cite  the  numerical  experiments  of  Canadas 
and  Kolb  (1986)  who  solved  this  problem  numerically  via  least-squares  data  fitting. 


FIGURE  CAPTIONS 


1.  Reference  and  perturbed  velocity  profiles  for  first  set  of  experiments  (Figures  1  -  5). 


2.  Reference  and  perturbed  wavelets  for  first  set  of  experiments. 


3.  Reference  and  perturbed  traces  for  p  =  0.0, 0.4. 


4.  Comparison  of  exact  and  recovered  velocity  profiles  from  data  of  Figure  3. 


5.  Comparison  of  exact  and  recovered  wavelets  from  data  of  Figure  3. 


6.  Perturbed  traces,  p  =0.0, 0.4;  1%  r.m.s.  perturbation  in  e(z),f(t).  Reference  traces  same  as 
in  Figure  3. 


7.  Comparison  of  exact  and  recovered  wavelets  from  data  of  Figure  6. 


8.  Comparison  of  exact  and  recovered  velocity  profiles  from  data  of  Figm;e  6. 


9.  Perturbed  traces,  p  =0.0, 0.4;  10 %r.m.s.  perturbation  in  c(z),f(t).  Reference  traces  same 
as  in  Figure  3. 


10.  Comparison  of  exact  and  recovered  wavelets  from  data  of  Figure  9. 


11.  Comparison  of  exact  and  recovered  velocity  profiles  from  data  of  Figure  9. 
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